
expr<-read.table("mRNA.txt",header = T,row.names = 1)
colnames(expr)<-gsub("\\.","-",colnames(expr))



group_list<-read.csv("sample_order_TCGA.csv",header = T,row.names = 1)

expr<-expr[,colnames(expr)%in%rownames(group_list)]
group<-data.frame(row.names = rownames(group_list),
                  group=group_list$group)

#group$group<-ifelse(group$group==1,"C1","other" )
#group$group<-ifelse(group$group==2,"C2","other" )
group$group<-ifelse(group$group==3,"C3","other" )
group<-as.matrix(group)

suppressMessages(library(limma))

design <- model.matrix(~0+factor(group))
colnames(design)=levels(factor(group))
rownames(design)=colnames(expr)
##两组以上的比较，寻找差异基因，使用topTableF函数
contrast.matrix<-makeContrasts("C3-other",
                               levels = design)

contrast.matrix ##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较


deg = function(expr,design,contrast.matrix){
  ##step1
  fit <- lmFit(expr,design)
  ##step2
  fit2 <- contrasts.fit(fit, contrast.matrix) 
  ##这一步很重要，大家可以自行看看效果
  
  fit2 <- eBayes(fit2)  ## default no trend !!!
  ##eBayes() with trend=TRUE
  ##step3
  tempOutput = topTable(fit2, coef=1, n=Inf)
  nrDEG = na.omit(tempOutput) 
  #write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
  head(nrDEG)
  return(nrDEG)
}

re = deg(expr,design,contrast.matrix)


final<-re
adj_pval_0.05_final<-final[which(final$P.Val<0.05),]
LogFC_1_final<-final[which(final$logFC>2),]
LogFC_f1_final<-final[which(final$logFC< -2 ),]
#取出上调基因的名字
up_gene<-intersect(row.names(LogFC_1_final),row.names(adj_pval_0.05_final))
down_gene<-intersect(row.names(LogFC_f1_final),row.names(adj_pval_0.05_final))
up<-length(up_gene)
down<-length(down_gene)
up_gene_file<-final[up_gene,]
down_gene_file<-final[down_gene,]
DEG_file<-rbind(up_gene_file,down_gene_file)

C3<-DEG_file

write.csv(DEG_file,file="C3_mRNA.csv")

unique(unique(rownames(C3),rownames(C1)),rownames(C2))




